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Abstract 

The  nonlinear  problem  of  tracking  and  predicting  where  a  satellite  will  be  over  some  time 
can  be  difficult  with  the  recognition  of  modeling  error  and  ground  site  radar  tracking  errors. 
For  this  reason  it  is  important  to  have  an  accurate  modeling  program  with  the  fidelity  to  cor¬ 
rect  for  any  errors  in  orbital  motion  and  predict  the  most  accurate  positioning  at  some  future 
time.  The  Extended  Kalman  Filter  is  one  such  program  that  can  accurately  determine  posi¬ 
tion  over  time  given  estimate  ranges  for  sources  of  error.  However,  the  Extended  Kalman 
Filter  contains  many  linear  approximations  that  allow  its  prediction  and  correction  methods 
to  work.  This  paper  will  discuss  the  effects  of  replacing  the  linearizing  approaches  made  in 
the  orbital  model  part  of  the  program  with  numerical  small-step  approaches.  The  overall 
errors  during  prediction  will  be  compared  for  an  analysis  of  the  corrective  ability  of  the 
filter.  Additionally  a  final  prediction  at  a  later  date  and  another  location  will  serve  as  an 
indicator  to  the  usefulness  of  the  prediction  capabilities  over  time. 

In  exploring  these  effects,  it  will  be  shown  that  the  linearizing  approximations  made  in 
the  development  are  a  good  approximation  to  the  numerical  results.  The  effects  of  modeling 
error,  perturbation  effects  included,  and  the  degree  of  approximation  all  play  a  significant 
role  in  accuracies  of  prediction.  The  effects  of  removing  linearizations  are  small  in  com¬ 
parison  to  the  effects  of  perturbations  and  modeling  error.  The  results  of  the  numerical  ap¬ 
proximations  contain  a  great  robustness  and  as  such  help  simplify  the  modeling  process.  The 
modeling  process  is  discussed  with  reference  to  ADA  program  code.  With  these  results  it 
can  be  seen  that  there  are  several  methods  of  using  the  Extended  Kalman  Filter  for  orbital 
prediction  which  maintain  a  high  degree  of  accuracy  and  can  be  very  useful  when  applied 
to  real-world  satellite  prediction. 

Introduction 

The  application  of  satellites  takes  many  forms  including  international  communi¬ 
cations  and  television  relay  to  space-based  telescopes.  However,  a  satellite’s  effec¬ 
tiveness  is  directly  based  on  the  ability  to  track  and  communicate  with  the 
spacecraft  over  its  lifetime.  In  order  to  track  a  satellite,  a  ground  station  will  use  a 
site-defined  coordinate  frame  to  track  the  motion.  This  coordinate  frame  uses  the 
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distance  from  the  site  to  the  satellite  (p),  the  azimuth  angle  (Az),  and  the  elevation 
angle  (El)  to  determine  where  in  the  sky  a  satellite  is  at  some  time.  The  ability  for 
a  ground  station  site  to  find  and  track  a  satellite  is  based  of  the  predictions  made  by 
mathematical  models  of  the  orbit.  The  more  accurate  a  model  can  be  at  predicting 
a  satellite’s  motion,  the  more  successful  a  ground  station  will  be  at  finding  and 
tracking  the  satellite. 

This  paper  discusses  methods  of  improving  an  Extended  Kalman  Filter  to  im¬ 
prove  the  orbital  determination  and  prediction  process.  Small  perturbations  felt  by 
an  orbiting  body  cannot  be  fully  modeled.  Additionally  errors  in  the  method  of  ob¬ 
taining  the  actual  position  and  velocity  from  radar  data  will  cause  errors.  For  these 
reasons  it  is  essential  to  use  a  filtering  device  such  as  the  Kalman  Filter  to  statisti¬ 
cally  determine  the  most  probable  position  of  a  satellite.  This  method  also  contains 
some  errors  in  its  prediction  algorithms  that  take  the  form  of  mathematical  lin¬ 
earizations.  By  seeking  methods  of  reducing  these  linearizations  this  paper  explores 
methods  of  improving  the  Extended  Kalman  Filter  for  better  overall  prediction  and 
orbit  determination.  This  approach  has  not  been  used  before  in  the  Extended 
Kalman  Filter. 

Because  there  are  only  a  limited  number  of  ground  stations  around  the  Earth,  a 
satellite  cannot  be  continuously  tracked.  Once  a  connection  is  lost,  only  the  pre¬ 
diction  of  where  the  satellite  will  be  at  some  later  time  can  help  in  reestablishing 
communication.  One  of  the  most  accurate  methods  of  establishing  a  model  that 
takes  the  orbital  determination  errors  into  account  is  the  Extended  Kalman  Filter. 
The  Extended  Kalman  Filter  (EKF)  is  a  stochastic  estimation  algorithm.  The  EKF 
simply  uses  a  weighted  statistical  average  of  the  difference  in  position  and  velocity 
inputs  predicted  from  the  model  and  known  from  the  ground  site  to  correct  the 
model  towards  more  precise  predictions  based  on  the  known  errors  in  those  inputs. 
The  EKF  can  be  tuned  to  use  anything  from  the  basic  two-body  orbit  model  to 
highly  accurate  multi-perturbation  model  for  its  predictions.  The  effects  of  recog¬ 
nizing  initial  modeling  error  can  also  be  explored  making  the  EKF  a  valuable  orbit 
prediction  and  modeling  tool.  Figure  1  shows  the  big  picture  of  how  the  EKF  is 
used  in  orbital  modeling.  Observations  from  one  site  allow  for  continual  orbital 
model  refinement  and  prediction  to  another  site. 


Extended  Kalman  Filter 
Big  Picture  Idea 


FIG.  1.  Big  Picture  Model  of  Orbital  Determination. 
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The  normal  Kalman  Filter  was  created  for  linear  applications.  Because  many 
systems  like  orbital  motion  are  nonlinear,  the  extended  filter  makes  many  approxi¬ 
mations  to  reduce  nonlinear  systems  to  linear  models.  For  these  reasons,  the  EKF 
has  been  the  subject  of  several  academic  studies.  At  the  Air  Force  Academy,  an 
EKF  was  designed  for  an  orbital  determination  scenario  and  the  filter’s  stability 
was  explored  over  a  wide  range  of  initial  inputs  and  models  [1].  The  same  initial 
EKF  orbit  model  and  scenario  is  explored  in  this  paper.  While  the  Academy’s 
effort  focused  on  the  stability  of  the  EKF  on  the  initial  orbit  model,  this  paper 
expands  from  this  to  explore  possible  methods  of  optimizing  the  EKF  results  by 
removing  linearizing  assumptions  made  during  the  initial  nonlinear  reduction  to  a 
linear  system. 

By  examining  two  methods  of  removing  linearizations,  this  paper  explores  the 
effects  of  linear  modeling  assumptions  as  well  as  what  inputs  have  the  greatest 
effects  on  the  orbit  prediction.  Linearizing  approximations  made  during  the  crea¬ 
tion  of  the  EKF  algorithm  are  replaced  with  numerical  approximations  that  allow 
for  nonlinear  perturbation  effects.  The  EKF  algorithm  attempts  to  minimize  the  dif¬ 
ference  in  predicted  and  observed  values  at  each  ground  site  observation.  To  ex¬ 
plore  the  effectiveness  of  removing  linearization  in  the  algorithm,  the  root  mean 
square  error  in  each  site  observation  in  position  (R)  and  velocity  (V)  is  explored 
over  a  pass  of  the  Mahe  Island  ground  station.  Fifteen  observations  every 
twenty  seconds  are  made  during  that  pass.  However,  the  true  prediction  power  of 
the  EKF  comes  from  not  only  accurately  predicting  the  next  observation  at  the 
same  ground  station,  but  also  accurately  predicting  an  observation  at  another  sta¬ 
tion  over  some  large  time.  The  errors  in  predicted  p,  Az,  and  El  at  the  Thule  ground 
station  are  also  explored  to  determine  the  effects  of  removing  linearizations.  The 
idea  of  model  refinement  towards  better  predictions  of  the  true  orbit  is  similar  to 
the  theoretical  mathematical  model  shown  later  in  Fig.  2.  The  effects  of  initial 
model  errors  and  the  addition  of  perturbations  to  the  model  are  also  explored  to  de¬ 
termine  the  effectiveness  of  removing  linearizations.  With  a  better  understanding  of 
the  issues  involved  in  orbital  determination  and  possible  ways  to  improve  these  re¬ 
sults,  conclusions  can  be  developed  on  the  linearizing  nature  of  the  Extended 
Kalman  Filter.  The  conclusions  developed  in  this  paper  will  allow  more  effective 
modeling  of  satellite  motion  for  the  insured  use  of  these  assets. 

The  Extended  Kalman  Filter 

Initial  Model  Assumptions 

The  Extended  Kalman  Filter  can  be  a  very  accurate  model  for  orbital  prediction; 
however,  there  are  several  assumptions  that  go  into  the  design  of  the  filter.  The  al¬ 
gorithm  itself  uses  a  basic  model  for  satellite  motion  based  on  the  two-body  equa¬ 
tion  of  motion.  Several  perturbation  effects  are  added  to  this  model  in  this  paper, 
even  though  there  is  no  way  to  fully  capture  all  of  the  characteristics  of  a  satellite’s 
motion.  The  very  nature  of  the  statistical  averaging  in  the  algorithm  also  makes  the 
assumption  that  none  of  the  input  data  into  the  system  are  completely  accurate  once 
any  error  is  specified  in  the  input.  As  the  filter  attempts  to  further  linearize  the  mo¬ 
tion  of  a  satellite  over  time,  a  first  or  second  order  Taylor  series  is  used  [2].  This  ap¬ 
proximation  makes  the  assumption  that  only  one  or  two  terms  of  a  Taylor  series  are 
a  good  approximation.  These  assumptions  help  to  linearize  the  nonlinear  problem 
of  orbital  motion.  The  purpose  of  this  paper  is  to  remove  some  of  the  linearization 
modeling  assumptions  in  the  Extended  Kalman  Filter. 
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Generalized  Math  Technique 

The  Extended  Kalman  Filter  algorithm  is  made  of  two  primary  pieces,  a  pre¬ 
diction  component  and  a  correction  component.  Both  of  these  components  help  to 
create  an  active  algorithm  that  will  statistically  compensate  a  model  of  an  orbit 
to  minimize  the  observed  and  predicted  position  values  and  velocity  state  values. 
The  EKF  requires  several  initial  inputs  in  order  to  work  properly.  These  inputs  in¬ 
clude:  the  latitude,  longitude,  and  altitude  of  the  ground  site;  the  individual  radar 
site  biases;  the  initial  error  covariance  matrix  (matrix  Pi);  the  error  expected  from 
the  radar  site  data  (vector  Rerror);  the  estimated  system  dynamic  modeling  error 
(matrix  Qt);  and  the  number  and  type  of  perturbation  effects  to  consider  [3].  The 
following  perturbations  are  explored:  J2,  J3,  J4,  drag,  Sun,  and  Moon  (the  satellite’s 
ballistic  coefficient  (BC)  will  be  given  to  calculate  the  drag  perturbation).  In  this 
project,  either  all  or  none  of  the  possible  perturbation  effects  are  considered.  The 
initial  Pi  matrix,  Rerror  vector,  radar  biases,  and  ground  site  data  are  given  for  this 
problem  and  will  remain  constant.  The  initial  Qt  is  either  zero  or  1  X  10”7.  (This  is 
the  same  modeling  error  estimation  assumed  from  the  original  Lyapunov  Stability 
analysis.)  The  nature  of  these  inputs  are  further  explored  throughout  this  paper.  The 
Qu  Rerror,  and  P,  quantities  are  shown  below  in  equations  (l)-(3).  The  six-by-six  na¬ 
ture  allows  for  the  position  and  velocity  in  each  inertial  axis  to  be  represented. 
Additionally,  the  initial  two-body  state  model  is  shown  in  equations  (4)  and  (5). 
With  these  inputs  and  the  initial  model  established,  the  EKF  is  posed  to  begin  its 
stochastic  estimation  process. 

Ri  0  0  0  0  0 

0  Rj  0  0  0  0 

0  0  Rk  0  0  0 

0  0  0  Vi  0  0 

0  0  0  0  Vj  0 

0  0  0  0  0  Vk 


Rerror  Crj  (2) 

crk 

Ri  0  0  0  0  0 

0  Rj  0  0  0  0 

0  0  Rk  0  0  0 

0  0  0  Vi  0  0 

0  0  0  0  Vj  0 

0  0  0  0  0  Vk 

X  =  R  (4) 

X=  V=  -4r+ Ap 
R2 


(5) 
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Prediction  Process 


The  Extended  Kalman  Filter’s  primary  purpose  is  to  establish  very  accurate  pre¬ 
dictions  of  an  orbit  over  time  knowing  the  error  that  is  expected  in  the  model  and 
the  actual  measurements  [1,2].  The  satellite’s  motion  is  initially  linearized  to  either 
a  first  or  second  order  Taylor  series  transition  matrix.  This  state  transition  matrix 
allows  for  the  prediction  of  the  states  of  position  (R)  and  velocity  (V)  at  some  fu¬ 
ture  time  (Ar).  The  state  transition  matrix  is  based  on  the  rate  of  change  in  V  verses 
the  rate  of  change  of  the  R  and  is  defined  as  the  F  matrix.  Both  the  state  transition 
matrix  (0)  and  the  F  matrix  are  shown  in  equations  (6)  and  (7)  respectively  [2]. 

A  t2 

=  e-m  =  /  +  f\t  +  F2—  (6) 


F 


ax 

ax 


(7) 


The  Cowell  method  of  orbit  propagation  is  used  to  update  the  states  and  the  error 
covariance  matrix  ( P ).  The  Cowell  method  of  orbit  determination  uses  the  fourth 
order  Runge-Kutta  (RK4)  approximation  method  to  solve  the  differential  equation 
for  the  new  states  (X)  at  some  time  (; tn ).  This  differential  equation  can  be  seen  below 
in  equation  (8),  where  variables  with  bars  over  the  top  are  the  predicted  values  and 
variables  with  hats  are  the  corrected  estimations.  The  RK4  method  propagates  the 
states  over  100  time  steps  between  every  20-second  observation  gap  and  between 
the  long  gap  in  time  to  the  Thule  site.  Error  covariance  of  the  new  states  will  be  pre¬ 
dicted  from  the  state  transition  and  the  estimated  error  covariance  as  shown  in  equa¬ 
tion  (9).  This  idea  can  be  seen  in  Fig.  2  showing  the  mathematical  process  of 
predicting  and  correcting  the  states  end  error  over  time  [2]. 


x(f„)  =  f  X  dr  +  mH) 

(8) 

P(t„)  =  (f)(tn,  t„-l)P(t„-l)<f)T(t„,  tn- 1) 

(9) 

From  estimated  values,  new  predicted  orbital  values  can  be  found,  but  these  val¬ 
ues  are  of  little  use  for  further  propagation  if  they  do  not  match  the  actual  observa¬ 
tions.  This  is  where  the  correction  process  of  the  EKF  identifies  its  usefulness.  The 
first  step  in  correcting  the  system  is  to  establish  the  error  between  predicted  values 
and  actual  observations.  This  is  done  through  the  H  matrix  as  shown  below  in 
equation  (10)  [2].  Because  it  is  known  that  the  radar  site  data  contains  errors  as  well 


FIG.  2.  Correction  Process. 
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as  the  predictions,  a  statistical  analysis  process  will  be  needed  to  correct  the  pre¬ 
dictions  and  create  new  estimations  of  the  states  [4].  The  next  step  in  the  process  is 
the  development  of  the  Kalman  Gain  matrix  ( K ).  The  K  matrix  uses  the  statistical 
weighed  average  of  the  predicted  error  covariance  and  the  known  radar  site  errors 
to  establish  a  gain  that  will  seek  to  minimize  the  diagonals  of  the  estimated  error 
covariance  matrix.  The  K  matrix  can  be  seen  in  equation  (11)  [2].  With  this  data, 
estimations  of  the  states  and  the  error  covariance  can  be  created  from  a  correction 
of  the  original  predictions  and  the  actual  site  information.  These  correction  equa¬ 
tions  can  be  seen  below  in  equations  (12)  and  (13),  where  the  [Z  —  gX]  term  is  the 
difference  in  the  observed  states  versus  the  predicted  states  [2]. 


w  dg(X(/„))  A  Observations 

m„)  ~  ax 

(10) 

K(tn)  =  P(t„)H'[H  P(t„)HT  +  Renot]-' 

(11) 

X(#„)  =  X(f„)  +  K[ Z  -  gX(tn)] 

(12) 

P{tn)  =  [/  -  KH]P(tn) 

(13) 

Orbit  and  Measurement  Performance 

Orbit 

The  EKF  scenario  that  is  explored  in  this  project  is  that  of  a  satellite  with  the  fol¬ 
lowing  set  of  orbital  elements: 


Semimajor  Axis: 

6697.0  km 

Eccentricity: 

0.0177 

Inclination: 

96.74° 

Long,  of  the  Ascending  Node: 

94.68° 

Argument  of  Perigee: 

57.35° 

Mean  Anomaly: 

0.0° 

The  given  satellite  is  the  same  as  that  used  in  the  original  Academy  study  [1].  This 
satellite  is  in  a  nearly  circular  low  Earth  orbit.  A  visual  representation  of  this  satel¬ 
lite  and  the  monitoring  ground  stations  can  be  seen  in  Fig.  3. 

Measuring  Sites  and  Performance  Data 

In  order  for  an  Extended  Kalman  Filter  to  work  as  an  orbital  optimization  tool,  it 
is  necessary  to  have  ground  site  observations  of  that  orbit  which  can  be  used  to 
correct  the  orbital  model.  Both  the  Mahe  and  Thule  ground  stations  are  used  for 
this  project  and  analysis.  The  Mahe  ground  station  is  located  at  4.8°  south  and 
55.5°  east.  This  station  will  provide  the  primary  measurements  of  the  p,  Az,  and  El 
of  the  satellite.  These  measurements  can  be  easily  converted  to  R  and  V  vectors  for 
the  satellite  given  a  known  location  and  time  of  the  observations.  The  date  of  these 
observations  is  11  May  1995.  The  biases  in  the  ground  site  radar  and  the  noise 
expected  from  these  observations  can  be  seen  in  Table  1 .  Computer  generated  obser¬ 
vations  of  the  orbit  with  error  can  be  seen  in  Table  2.  The  RMS  error  of  these 
observations  to  the  predicted  results  at  each  observation  is  used  as  one  measure  of 
performance.  The  orbital  prediction  capability  of  the  EKF  can  then  be  compared  to 
another  ground  site  to  determine  the  accuracy  of  the  model.  The  second  site  used 
for  observation  is  that  of  Thule,  Greenland  at  76.57°  north  and  68.3°  west.  The 
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FIG.  3.  STK  Graphical  Representation  of  Modeling  Scenario. 


biases  and  expected  noise  of  this  site  can  be  seen  in  Table  3.  This  observation  is 
made  21  minutes  and  20  seconds  after  the  last  Mahe  observation.  Error  difference 
between  this  observation  and  model’s  predicted  results  will  be  compared  to  provide 
another  measure  of  performance  [4,  5]. 


TABLE  1.  Measurement  Characteristics  of  Mahe  Ground  Station 


Measurement 

Biases 

Standard  Deviation  a 

P  (km) 

0.15 

0.15 

Az  (deg) 

0.0001 

0.01 

El  (deg) 

0.0001 

0.01 

TABLE  2. 

Measurement  Data  from  Mahe  Ground  Station 

Hour 

Minutes 

Seconds 

r  (km) 

Az  (deg) 

El (deg) 

11 

50 

0 

1770.334 

195.1388 

1.737 

11 

50 

20 

1636.96 

197.6087 

3.0125 

11 

50 

40 

1507.572 

200.5082 

4.414 

11 

51 

0 

1382.047 

203.9435 

5.8441 

11 

51 

20 

1262.499 

208.0942 

7.4128 

11 

51 

40 

1150.56 

213.1119 

9.0443 

11 

52 

0 

1049.464 

219.2097 

10.7308 

11 

52 

20 

961.3957 

226.5796 

12.3925 

11 

52 

40 

891.6892 

235.4314 

13.8777 

11 

53 

0 

844.4838 

245.5759 

14.9465 

11 

53 

20 

823.1501 

256.6734 

15.3395 

11 

53 

40 

829.9796 

267.9568 

14.9301 

11 

54 

0 

865.0305 

278.6099 

13.8245 

11 

54 

20 

924.6791 

288.0676 

12.2013 

11 

54 

40 

1003.953 

296.0857 

10.3958 
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TABLE  3.  Measurement  Characteristics  of  Thule  Ground  Station 


Measurement 

Biases 

Standard  Deviation  <x 

P  (km) 

0.0708 

0.026 

Az  (deg) 

0.0013 

0.026 

El  (deg) 

0.0075 

0.022 

Matching  Results 

In  the  analysis  of  an  EKF,  it  is  important  to  establish  an  accurate  reference  model 
from  which  to  expand.  First,  the  results  of  the  initial  Academy  study  were  dupli¬ 
cated  using  the  EKF  in  [1].  Several  more  accurate  constants  of  motion  and  pertur¬ 
bation  effects  were  added  to  this  model.  The  results  of  the  Academy’s  study  were 
reproduced  and  a  comparison  of  these  results  can  be  seen  in  Table  4.  Matching 
these  results  helped  establish  a  baseline  for  improvement. 

Nominal  Results 

The  following  figures  represent  the  possible  input  options  and  the  error  results 
for  the  Mahe  observations.  Some  of  the  possible  input  options  on  the  initial  EKF 
are  as  follows:  1st  versus  2nd  order  Taylor  series,  zero  versus  all  perturbation  ef¬ 
fects,  and  zero  versus  1  X  10”7  modeling  error  for  the  filter.  It  can  be  seen  from 
these  results  that  the  largest  difference  in  model  results  can  be  seen  with  the  addi¬ 
tion  of  modeling  error.  Both  Figs.  4  and  5  show  these  results  for  the  RMS  position 
and  velocity  results  respectively.  It  is  important  to  note  that  all  of  the  RMS  error 
plots  in  this  report  are  graphed  with  respect  to  the  Julian  date  (JD)  of  observations. 
However,  these  observations  are  made  twenty  seconds  apart  as  seen  in  Table  2. 

Design  Improvement  Methods 

The  focus  of  this  project  is  on  the  improvement  of  the  EKF  orbital  determination 
modeling  process.  In  particular,  improving  the  EKF  by  removing  linearizing  as¬ 
sumptions  will  be  explored.  Assumptions  on  F  and  in  particular  are  the  interest 
of  this  analysis.  These  matrices  are  concerned  with  the  prediction  portion  of  the 
EKF  [6,  7].  Improving  the  prediction  process  is  one  method  of  improving  the  over¬ 
all  modeling  capability.  There  are  many  other  methods  of  improvement  that  can  be 
explored,  however  the  biggest  advantage  in  modeling  capability  would  be  the  ad¬ 
dition  of  all  the  known  perturbation  effects  and  the  removal  of  these  linearization 
results.  For  each  of  the  linearizing  assumptions,  the  mathematical  techniques  used 
to  remove  the  assumption  and  the  results  of  these  changes  are  discussed  in  follow¬ 
ing  sections.  Many  of  the  outlining  causes  of  these  results  will  be  analyzed.  Con¬ 
clusions  on  possible  uses  of  these  results  and  nonlinearizations  are  discussed. 


TABLE  4.  EKF  Data  Matching  Results  Comparison:  At  Thule’s  First  Observation  Using 
Original  1995  Data 


Observed 

Computed 

Calculated  Error 

Error  Results 

Rho  (km) 

1194.8692 

1205.322457 

10.453257 

38.54 

Az  (degs) 

34.8631 

35.53747 

0.67437 

1.38 

El  (degs) 

6.3154 

6.569847 

0.254447 

0.78 
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0.24 


0.22  - 


0.2 


0.18 


0.16 


0.14 


0.12 


0.1 
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* 


Both  Sets  of  Phi  =  1 


2.4498  2.4498  2.4498  2.4498  2.4498  2.4498 
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FIG.  4.  Initial  RMS  Error  Plot  of  Position  (km)  Including  Modeling  Error. 


JD  Time  (20  sec  steps) 

FIG.  5.  Initial  RMS  Error  Plot  of  Velocity  (km/sec)  Including  Modeling  Error. 
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F  Matrix  Nonlinearization 

In  the  initial  development  of  Extended  Kalman  Filter  model,  the  change  in  the 
rates  of  change  in  the  states  versus  the  change  in  states  themselves,  or  F  matrix,  is 
important  to  the  understanding  and  modeling  of  orbital  motion.  The  initial  F  ma¬ 
trix  equation  can  be  seen  in  equation  (7).  This  equation  is  a  very  key  component  to 
developing  orbital  motion  and  the  EKF  algorithm.  However,  in  order  for  the  EKF 
to  use  this  matrix,  an  analytic  approximation  approach  must  be  used  to  reduce  the 
equations  of  motion  into  a  usable  form.  As  seen  in  equation  (14),  most  parts  of  the 
matrix  go  to  zero  or  to  identity  except  for  one  portion.  Using  only  the  two-body 
equations  motion,  analytic  approximations  for  each  of  the  terms  in  nontrivial  por¬ 
tion  of  equation  (14)  can  be  found  [2],  [3],  [5].  This  lower  nontrivial  portion  can  be 
seen  in  equation  (15).  A  series  of  nine  equations  can  be  painstakingly  developed  to 
solve  for  the  final  F  matrix.  However  in  the  development  of  this  approach  to  lin¬ 
earizing  the  F  matrix  differentials,  many  of  the  effective  perturbation  terms  are  lost 
and  cannot  be  used  in  the  orbital  motion  model  [2]. 

ait  1 


av 

aR 

The  two-body  acceleration  is  denoted  by  the  components  Ai,  Aj,  and  Ak. 
Mathematical  Technique  and  Implementation 


In  order  to  remove  this  linearization  effect  and  include  all  of  the  perturbation 
terms  in  the  model  another  approach  to  creating  the  F  matrix  must  be  developed. 
The  equation  for  this  development  can  be  seen  in  equation  (16)  below.  The  actual 
implementation  of  this  equation  requires  a  change  in  the  position  vector  that  is  then 
feeding  the  change  back  into  the  RK4  prediction  algorithm  to  determine  the  accel¬ 
erations  that  would  have  been  felt  by  that  change.  Small  changes  in  each  inertial  di¬ 
rection  (i,j,  and  k)  for  each  term  will  complete  the  lower  portion  of  equation  (15). 
By  assuming  only  small  perturbations  in  the  model,  the  small  numerical  changes  in 
the  rate  of  change  in  states  can  be  compared  to  the  small  numerical  changes  in  the 
states  themselves.  This  approximation  is  very  similar  to  the  finite  differential  analy¬ 


sis  development  of  the  H  matrix. 


I 

0 


Where 


AV 

AR 


A Ai 

A Ai 

A Ai 

A  Ri 

A  Rj 

ARk 

A Aj 

A Aj 

A Aj 

A  Ri 

ARj 

ARk 

A Ak 

A Ak 

A Ak 

A  Ri 

ARj 

ARk 

(16) 
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Note  that  the  acceleration  now  includes  two-body  accelerations  as  well  as  all  the  re¬ 
quested  perturbation  accelerations. 

In  the  development  of  this  approach  it  is  necessary  to  understand  what  sort  of 
“small”  step  in  position  will  accurately  model  the  acceleration  perturbations  that 
may  be  felt.  To  get  an  idea  of  the  needed  step  size  in  position  that  should  be  taken, 
the  acceleration  felt  by  the  system  was  graphed  with  respect  to  the  position.  These 
results  can  be  seen  in  Figs.  6  and  7.  By  inspection,  a  position  step  size  of  0.001  R, 
roughly  6  km,  should  be  able  to  accurately  model  the  perturbations.  This  is  because 
a  small  size  is  needed  to  model  the  curl  at  the  end  of  the  acceleration  in  the  j  di¬ 
rection.  However,  to  ensure  these  results,  the  step  size  was  tested  over  a  wide  range 
as  seen  in  the  comparison  figures. 

Results 

The  result  of  this  implementation  in  the  Extended  Kalman  Filter  algorithm  can 
be  seen  in  the  sets  of  data  from  both  the  Mahe  and  Thule  ground  stations.  The  RMS 
error  in  both  R  and  V  were  plotted  versus  time  of  the  observation.  These  results 
compared  with  this  nominal  model  results  and  the  same  inputs  can  be  seen  in  the 
figures.  As  can  be  seen  from  these  results,  there  is  relatively  no  difference  in  the 
modeling  ability  of  the  new  numerical  F  approach.  These  results  are  displayed  in 
Figs.  8  and  9.  To  explore  what  effects  differing  step  sizes  would  have  a  range  of 
multiplication  factors  were  input  into  the  same  model  as  can  be  seen  in  Fig.  10.  Due 
to  the  fact  that  the  small  position  changes  will  only  affect  velocity  terms,  only  the 
RMS  velocity  errors  will  be  shown  to  explore  different  step  sizes.  These  results  can 
be  seen  in  Fig.  10.  Figure  1 1  goes  on  to  explore  the  effects  of  modeling  error  on  the 
numerical  F  matrix,  again  only  the  velocity  is  shown  because  that  is  the  only 
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FIG.  6.  Acceleration  Versus  Position  Curve  in  the  i  and  k  Directions. 
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FIG.  7.  Acceleration  Versus  Position  Curve  in  the  j  Direction. 
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FIG.  8.  Numerical  F  Matrix  RMS  Error  Plot  Comparison  of  Position  (km). 
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FIG.  11.  Numerical  F  Matrix  RMS  Error  Plot  of  Velocity  with  Modeling  Error  (Step  Factor  =  0.001). 


factor  that  experiences  any  observable  effect.  It  can  also  be  seen  from  these  graphs 
that  using  a  second  order  Taylor  series  approximation  will  yield  results  similar  to 
the  initial  results,  which  include  modeling  error.  However,  using  only  a  first  order 
Taylor  series  will  negate  the  effects  of  any  modeling  error.  The  best  results  come 
from  the  largest  step  size  multiplication  factor,  roughly  60  km.  These  results  are 
very  small  and  have  virtually  no  effect  on  the  system.  For  all  purposes  the  step  size 
is  not  that  significant  of  a  factor  in  the  F  matrix  approximation. 

To  examine  the  effects  of  this  F  matrix  correction,  several  inputs  were  varied  and 
the  RMS  outputs  were  plotted.  These  results  from  varying  the  Taylor  series  ap¬ 
proximation  are  used  and  the  number  of  perturbation  effects  can  be  seen  in  the  pre¬ 
vious  Figs.  10  and  11.  The  results  appear  to  be  very  similar  to  those  of  the  previous 
cases.  However  there  are  larger  and  more  variable  RMS  error  results  with  the  second- 
order  Taylor  series.  Upon  inspection  of  all  of  the  final  Thule  results  of  the  test  cases 
with  their  variable  changes,  an  interesting  fact  arises.  The  cases  where  all  of  the 
perturbation  effects  were  added  to  the  model  produce  better  results.  However,  the 
cases  that  include  a  second-order  modeling  term  are  far  more  accurate  than  all  of 
the  above.  These  error  results  in  the  predicted  versus  actual  p,  Az,  and  El  can  be 
seen  in  Table  5.  This  table  includes  all  of  the  initial  test  cases  with  the  original  EKF 
as  well  as  all  of  the  test  results  with  the  new  numerical  F  matrix  approach.  It  is  also 
important  to  note  that  there  was  no  effect  on  these  results  with  the  implementation 
of  modeling  error  when  a  first  order  Taylor  series  approximation  was  used. 

Numerical  F  Matrix  Conclusions 

It  can  be  seen  from  the  results  above  that  the  largest  variation  in  modeling  results 
comes  from  the  accuracy  of  the  order  of  the  </>  matrix  Taylor  series  approximation. 
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The  use  of  the  second  order  Taylor  series  is  accurate  enough  to  actually  capture  the 
effects  of  a  modeling  error.  The  results  from  each  test  case  nearly  match  all  of  the 
initial  data  results  with  similar  initial  configurations.  In  cases  where  more  error 
arose  than  the  initial  model,  it  appears  that  more  rounding  error  played  a  part  in 
looking  at  the  small  acceleration  changes.  The  F  matrix  still  relies  on  the  RK4  al¬ 
gorithm,  which  inherently  contains  some  error.  In  conclusion,  it  seems  that  the 
F  matrix  two-body  linearization  is  a  good  assumption  to  make,  and  the  largest  dif¬ 
ferences  come  in  the  manipulation  of  the  state  transition  matrix.  The  results  using 
a  better  </>  matrix  approximation  are  significantly  greater  and  thus  beg  for  a  better 
understanding.  For  this  reason  completely  removing  the  Taylor  series  approxima¬ 
tion  and  the  need  for  the  F  matrix  from  the  c/>  matrix  is  explored. 

cf)  Matrix  Nonlinearization 

The  matrix  is  the  state  transition  matrix  of  the  Extended  Kalman  Filter.  This 
matrix  allows  the  filter  to  predict  how  the  error  covariance  of  the  states  will  change 
with  time  [2].  In  the  current  use  of  the  EKF,  the  <fi  matrix  is  approximated  by  a  sec¬ 
ond  order  Taylor  series  as  seen  in  equation  (6).  However,  this  approximation  can 
cause  error  in  the  loss  of  the  higher  order  Taylor  series  terms  used  to  solve  the  dif¬ 
ferential  equation  shown  in  equation  (17).  This  equation  is  from  what  derives  the 
initial  need  for  the  F  matrix.  By  using  a  different  approach,  the  idea  of  a  state  tran¬ 
sition  matrix  can  be  solved  without  the  need  of  a  differential  equation,  and  thus 
without  the  need  of  the  Taylor  series  or  the  F  matrix. 

<t>  =  F(f>  (17) 

Mathematical  Technique  and  Implementation 

The  <fi  matrix  is  merely  the  equations  of  state  necessary  to  transfer  from  the  ini¬ 
tial  states  to  some  final  state  [6].  This  relationship  can  be  seen  in  equation  (18)  [2]. 
Much  like  the  numerical  finite  differential  analysis  approach  taken  with  F,  the  </> 
matrix  will  perturb  an  initial  state  and  propagate  it  forward  to  compare  it  with  the 
nominal  propagated  state.  By  establishing  the  difference  in  perturbations  results  in 
each  direction  compared  to  the  initial  perturbations,  the  </>  matrix  can  be  directly 
solved  for  as  seen  in  equation  (19). 

8X  =  (f)8X  (18) 
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Results 

The  results  of  this  simple  algorithm  are  very  elegant  and  also  very  surprising.  In 
the  establishment  of  the  new  algorithm,  it  was  once  again  necessary  to  determine  a 
proper  step  size  in  position.  The  ratio  of  the  position  step  size  to  the  velocity  step 
size  is  also  a  variable  to  be  determined.  The  ratio  value  of  0.001  (which  was  used 
in  previous  studies)  seemed  the  best  choice  since  the  velocity  was  on  average  three 
orders  of  magnitude  less  than  the  position.  The  effects  of  perturbations  in  the  orbit 
as  well  as  the  effects  of  modeling  error  were  explored  in  this  analysis.  The  com¬ 
parison  of  these  variables  with  the  initial  estimate  at  position  step  size  ( dx)  and  ratio 
factors  can  be  seen  in  Figs.  12  and  13.  The  greatest  effects  on  the  data  in  this  case 
occur  when  the  modeling  error  is  applied.  When  compared  to  the  results  from  the 
initial  filter  plots,  these  plots  are  nearly  identical  (this  result  can  also  be  seen  in  the 
Thule  data  site  comparisons).  The  use  of  modeling  error  matches  the  second-order 
Taylor  series  approximation  and  without  the  use  of  modeling  error  matches  the 
first-order  results.  These  initial  results  can  be  found  in  Figs.  4  and  5.  In  the  com¬ 
parison,  when  looking  at  the  RMS  error  velocity  plots  of  varying  dx  and  ratio  val¬ 
ues  over  several  orders  of  magnitude,  it  can  be  seen  that  there  is  hardly  any  change 
in  prediction  errors.  These  comparison  graphs  can  be  seen  in  Figs.  14  and  15  re¬ 
spectively.  Once  again  only  the  velocity  graphs  are  shown  here  because  they  show 
the  actual  correlation  to  the  new  changes  in  adding  acceleration  terms.  These  results 
show  that  the  new  numerical  </>  approach  is  very  robust  and  eliminates  the  need  for 
precise  refinement  when  it  comes  to  the  dx  and  ratio  estimations. 

The  final  prediction  data  produced  from  all  of  the  possible  test  cases  at  the  Thule 
ground  can  be  seen  in  Table  6.  This  data  also  includes  a  more  robust  set  of  varied 
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FIG.  12.  Numerical  </>  Matrix  RMS  Error  Plot  Comparison  of  Position  (km). 
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FIG.  13.  Numerical  </>  Matrix  RMS  Error  Plot  Comparison  of  Velocity  (km/sec). 
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FIG.  14.  Numerical  </>  Matrix  RMS  Error  Plot  of  Velocity  (km/sec)  with  dx  Step  Factor. 


An  Improved  Kalman  Filter  for  Satellite  Orbit  Predictions 


19 


5.2 


5.1 


x  10 


5  - 


Error 


4.9 


4.8  - 


4.7 


4.6 


4.5 


4.4 


*  * 

4 

. . : . 

* 

* 

*  Ratio  =  0.00001 

*  Ratio  =  0.001 

*  Ratio  =  100.0 


4.3 

2.4498  2.4498  2.4498  2.4498  2.4498  2.4498  2.4498  2.4498  2.4498 

,r.-r  ,  XlO6 

JD  Time  (20  Sec  Steps) 

FIG.  15.  Numerical  </>  Matrix  RMS  Error  Plot  of  Velocity  (km/sec)  with  Ratio  Step  Factor  Comparison. 


dx  and  ratio  results  and  can  be  compared  to  the  initial  EKF  results  found  in  Table  5. 
In  this  table,  it  can  be  seen  that  the  first  four  cases  involving  the  perturbations  and 
the  modeling  error  nearly  match  the  initial  model  results  with  the  same  variable 
changes.  These  results  help  to  show  that  the  method  works  in  its  approximation  and 
that  it  is  slightly  more  consistent  than  the  original  approximation.  The  unusual  re¬ 
sult  from  this  entire  test  is  in  the  fact  that  the  model  is  very  robust  at  handling  a 
large  range  in  position  and  velocities  step  sizes.  The  model  was  tested  to  breaking, 
failing  the  code,  in  the  maximum  (925  step  factor)  and  minimum  (1  X  10”7  step 
factor),  in  allowable  ranges  of  dx  and  in  ratio  step  size  multiplication  factors.  These 
tests  produced  no  significant  improvement  in  results  other  than  the  proof  that  the 
model  is  very  robust.  In  looking  at  the  direct  compensation  error  value,  [Z  —  gX], 
that  goes  into  the  Kalman  gain  matrix  to  correct  the  system  after  each  prediction, 
the  filter  worked  hard  to  correct  the  difference  down  to  zero.  Again  there  was  no  re¬ 
sponse  to  varying  the  dx  value  used.  These  results  are  shown  in  Fig.  16.  The  pre¬ 
dicted  mathematic  modeling  process  from  Fig.  2  is  reminiscent  of  the  actual 
process  shown  in  Fig.  16. 

Numerical  </>  Matrix  Conclusions 

It  can  be  seen  from  these  results  that  the  approach  of  removing  the  linearizing 
approximations  in  the  </>  matrix  and  allowing  for  all  of  the  perturbation  effects  al¬ 
lows  the  EKF  to  become  very  robust  at  modeling  data  of  any  sort.  While  it  is  still 
not  completely  understood  how  the  EKF  can  handle  such  a  wide  range  in  step  size 
inputs,  it  is  plain  to  see  that  the  best  results  come  from  the  use  of  modeling  error 
and  including  all  of  the  possible  perturbation  effects.  This  approximation  is  very 
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FIG.  16.  Numerical  f  Matrix  Direct  Compensation  Error  Value,  [Z  —  gX]. 


easy  to  implement  in  code  and  serves  as  a  great  modeling  source  as  it  is  consistently 
better  than  any  other  approximation  method.  This  <fi  matrix  fix  should  be  imple¬ 
mented  in  future  uses  to  improve  the  filter  operations.  This  is  because  it  is  con¬ 
stantly  accurate  while  allowing  for  a  robust  range  of  step  inputs  and  the  use  of  an 
imperfect  model. 

Overall  Conclusions  and  Recommendations 

It  has  always  been  the  pursuit  of  science  and  engineering  to  try  and  find  the  best 
way  to  produce  the  most  accurate  results.  In  the  case  of  the  Extended  Kalman  Fil¬ 
ter,  it  is  very  hard  to  create  the  better  mousetrap  of  orbital  determination  and  mod¬ 
eling.  The  Kalman  filter  is  not  just  a  simple  process  but  involves  a  variety  of 
external  inputs,  complex  matrix  operations,  and  statistical  determination  to  create 
an  optimizing  algorithm  for  orbital  determination.  In  the  approach  of  removing  lin¬ 
earizing  assumptions  made  during  the  initial  creation  of  the  EKF,  it  can  be  seen  that 
the  results  are  not  significantly  better  than  those  with  the  linear  approximations.  It 
appears  that  the  greatest  improvement  in  results  comes  from  knowing  how  to  set  up 
the  filter  and  what  inputs  are  necessary.  In  looking  at  the  comparison  results  of  each 
method,  it  appears  that  the  model  inputs  used  have  the  largest  effects  on  the  system. 
These  results  also  serve  to  prove  the  fact  that  the  initial  linearizations  and  approx¬ 
imations  created  are  a  fairly  accurate  representation  of  the  actual  process.  In  light 
of  the  results  of  the  </>  matrix  improvement,  it  is  recommended  that  this  change  be 
added  in  future  systems.  This  matrix  allows  for  a  decrease  in  needed  inputs  into  the 
model  and  allows  for  a  wide  range  of  possible  step  multiplication  factors.  This  sys¬ 
tem  would  help  ease  implementation  into  orbital  determination  problems  where  the 
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modeling  accuracy  is  never  truly  zero.  While  the  perfect  orbital  determination  soft¬ 
ware  can  never  be  developed  to  handle  all  satellites,  the  Extended  Kalman  Filter 
does  a  remarkable  job  of  orbital  modeling  given  some  initial  observations  for  fur¬ 
ther  optimization. 

This  approach  to  removing  some  of  the  linearization  in  the  system  model  is  not 
unique  to  the  orbital  determination  problem.  It  is  general  enough  that  it  can  be  ap¬ 
plied  to  the  estimation  of  almost  any  nonlinear  system. 
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